5 de Mayo de 2018
Entender cómo usar los datos geolocalizados para problemas de marketing
En especial
Explicar la demografía de los visitantes de cada centro comercial
Usar los datos de movimientos de la gente para caracterizar la atracción de cada centro comercial
Estudiar problemas de logística/cierre centros comerciales
Recordamos: geomarketing es un tipo de "location intelligence" que, utilizando datos geolocalizados de tiendas, clientes y competencia permite:
El objetivo de este Taller es encontrar las zonas de atracción en los barrios de Madrid para cada uno de los Ikeas (u otros centros comerciales)
Para ello:
Vamos a utilizar datos geolocalizados de Twitter para determinar donde viven los usuarios que van a cada uno de los centros de Ikea de Madrid
Hacemos un modelo (fit) para determinar mediante el modelo de Huff cual es la probabilidad \(p_{ij}\) de que un usuario en un barrio \(i\) visite el centro de Ikea \(j\)
Con ese modelo determinamos las zonas de atracción de cada Ikea, asignando cada barrio \(i\) a un Ikea \(j\) si la probabilidad \(p_{ij}\) es mayor de las \(p_{i*}\)
Vamos a utilizar una base de datos de tweets geolocalizados durante 2 años
load("./data/tweets_madrid.RData")
head(tweets_madrid)
## lat lon user time ## 21 40.49162 -3.689221 1 2012-12-01 00:03:40 ## 22 40.41296 -3.723460 2 2012-12-01 00:05:52 ## 36 40.41916 -3.706140 3 2012-12-01 00:03:07 ## 48 40.39049 -3.725962 4 2012-12-01 00:04:01 ## 54 40.42712 -3.673200 5 2012-12-01 00:07:19 ## 63 40.40596 -3.683690 6 2012-12-01 00:04:22
Cargamos los shapefiles de los barrios de Madrid
map_barrios <- readOGR(dsn ="./data/Barrios Madrid/",layer="200001469",
encoding = "latin1",stringsAsFactors=F)
## OGR data source with driver: ESRI Shapefile ## Source: "./data/Barrios Madrid/", layer: "200001469" ## with 128 features ## It has 3 fields
Los transformamos a la proyección WSG84
WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
mapb <- spTransform(map_barrios,CRSobj=WGS84)
proj4string(mapb)
## [1] "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
Los visualizamos
leaflet(mapb) %>% addTiles() %>% addPolygons(weight=1)
También necesitaremos la población de los barrios de Madrid
http://www-2.munimadrid.es/TSE6/control/seleccionDatosBarrio
pob_barrios <- read.csv("./data/Madrid-Poblacion_Barrios.csv")
Añadimos a la tabla del mapa una columna que sea la población (uniendo por el GEOCODIGO)
mapb@data$pob <- pob_barrios$Población[match(as.numeric(mapb$GEOCODIGO),
pob_barrios$GEOCODIGO)]
head(mapb@data)
## CODBDT GEOCODIGO DESBDT pob ## 0 797664 079011 011 Palacio 22471 ## 1 797665 079012 012 Embajadores 46264 ## 2 797666 079013 013 Cortes 10717 ## 3 797667 079014 014 Justicia 16847 ## 4 797668 079015 015 Universidad 31714 ## 5 797669 079016 016 Sol 7577
En concreto los centros comerciales en la comunidad de Madrid (http://www.madrid.org/nomecalles/)
centros <- read.table("./data/centros_comerciales.csv",header=T,stringsAsFactors = F)
head(centros,2)
## CMUN ETIQUETA MUNICIPIO ## 1 5 Leroy Merlin Alcalá de Henares ## 2 7 Leroy Merlin Alcorcón ## DIRECCIÓN ## 1 Ctra Nacional II, km. 34 - Calle Tamajón s/n (Centro Comercial La Dehesa) ## 2 Avenida de Europa, 26 (Centro Comercial Parque Oeste - Ctra. Nacional V, km. 14) ## lon lat ## 1 -3.327188 40.50292 ## 2 -3.849481 40.34330
Estimamos donde está el domicilio de cada usuario a nivel de barrio: barrio más frecuentado para hacer los tweets.
Para ello añadimos una columna GEOCODICO en la tabla de tweets de cada barrio:
pto_matrix <- cbind(tweets_madrid$lon,tweets_madrid$lat) pps1 <- SpatialPoints(pto_matrix,proj4str=WGS84) oo1 <- over(pps1,as(mapb,"SpatialPolygons")) tweets_madrid$GEOCODIGO <- mapb@data$GEOCODIGO[oo1] head(tweets_madrid)
## lat lon user time GEOCODIGO ## 21 40.49162 -3.689221 1 2012-12-01 00:03:40 079086 ## 22 40.41296 -3.723460 2 2012-12-01 00:05:52 079102 ## 36 40.41916 -3.706140 3 2012-12-01 00:03:07 079016 ## 48 40.39049 -3.725962 4 2012-12-01 00:04:01 079112 ## 54 40.42712 -3.673200 5 2012-12-01 00:07:19 079042 ## 63 40.40596 -3.683690 6 2012-12-01 00:04:22 079031
Nos quedamos solo con los tweets que suceden en alguno de los barrios
tws <- tweets_madrid[!is.na(tweets_madrid$GEOCODIGO),]
Es posible que haya usuarios que han venido a Madrid de vacaciones y hayan hecho muchos tweets. Nos quedamos solo con los tweets de aquellos usuarios que al menos haya estado en Madrid más de 5 meses. También nos quedamos sólo con aquellos que hayan hecho más de 10 tweets y menos de 10000.
tws$month <- month(tws$time)
users_table <- tws %>%
group_by(user) %>%
summarize(ntot=length(user),nmonths=length(unique(month)))
users_sel <- users_table$user[users_table$nmonths > 5 & users_table$ntot > 10 & users_table$ntot < 10000]
tws <- tws[tws$user %in% users_sel,]
Finalmente calculamos cual ha sido el barrio más frecuentado y cuánto ha sido frecuentado. Para cada usuario obtenemos cuantos tweets tiene en cada barrio
user_barrio_tweets <- tws %>% group_by(user,GEOCODIGO) %>% summarize(nt = length(user))
Ordenamos y calculamos el barrio más visitado
user_barrio_tweets <- user_barrio_tweets[order(user_barrio_tweets$user,user_barrio_tweets$nt,decreasing=T),]
user_barrio_most <- user_barrio_tweets %>%
group_by(user) %>%
summarize(ntot = sum(nt),
nmax=head(nt,1),
bmax=head(GEOCODIGO,1),
nbarrios=length(unique(GEOCODIGO)))
Solo consideramos aquellos en los que más del 30% de los tweets suceden en el barrio más visitado
user_barrio_most <- filter(user_barrio_most,ntot < 10000 & ntot > 10 & nmax/ntot > 0.3 & nbarrios > 3)
Vamos a ver cómo se comparan el número de homes con la población de cada barrio
aa <- aggregate(user_barrio_most$bmax,by=list(user_barrio_most$bmax),length)
colnames(aa) <- c("GEOCODIGO","nusers")
aa$GEOCODIGO <- as.numeric(aa$GEOCODIGO)
mm <- merge(aa,pob_barrios)
plot(mm$nusers,mm$Población,xlab="Twitter users",ylab="Población")
Estimamos cuantos usuarios han ido a un centro comercial, desde que barrio y a qué distancia.
Vamos a probar con los de Ikea. Con el de Vallecas
centros_sel <- centros[centros$ETIQUETA=="Ikea",] lon_Ikea <- centros_sel[3,"lon"] lat_Ikea <- centros_sel[3,"lat"]
Vemos qué usuarios han geolocalizado un tweet a menos de 500 metros del centro de IKEA
distancia_a_Ikea <- spDistsN1(cbind(tws$lon,tws$lat),c(lon_Ikea,lat_Ikea),longlat = T) users_Ikea <- unique(tws$user[distancia_a_Ikea < 0.5])
Hacemos una tabla de los barrios de donde vienen esos usuarios y la distancia (al centro del barrio)
users_barrio_Ikea <- user_barrio_most[user_barrio_most$user %in% users_Ikea,] agg_barrio_Ikea <- aggregate(users_barrio_Ikea$bmax,by=list(users_barrio_Ikea$bmax),length)
Juntamos estos datos con el mapa
mapb_Ikea <- mapb mapb_Ikea@data$nusers <- agg_barrio_Ikea$x[match(mapb_Ikea@data$GEOCODIGO,agg_barrio_Ikea$Group.1)]
Pintamos un choropleth con estos datos
mapp <- mapb_Ikea[!is.na(mapb_Ikea@data$nusers),]
xx <- mapp@data$nusers
xx <- (xx-min(xx))/(max(xx)-min(xx))
m <- leaflet(data = mapp) %>% addTiles() %>%
addPolygons(fillColor= ~ rgb(xx,0,0),
weight=1,fillOpacity = 0.8,
popup=~paste(DESBDT,"<br> nusers = ",nusers,sep=""))
Veamos como es este choropleth
m
Ejercicio 1:
¿Desde qué el barrio van más personas al centro de Ikea del norte?
¿Cuál es la distancia media que recorren los clientes para ir a cada uno de los Ikeas?
Finalmente, vamos a modelizar la atracción a cada uno de los centros comerciales (Ikea). Para ello vamos a utilizar la ley de Huff que dice que la probabilidad de que una persona que vive en \(i\) visite un centro comercial en \(j\) viene dada por
\[ p_{ij} \sim \frac{A_j}{d_{ij}^\beta} \]
donde \(A_j\) es lo atractivo que es el centro comercial \(j\) (normalmente se toma el tamaño), mientras que \(d_{ij}\) es la distancia del cliente \(i\) al centro \(j\).
En nuestro caso vamos a empezar por ver cómo la distancia determina el número de usuarios que va a un centro comercial, mediante la ley
\[N_{ij} \sim \frac{1}{d_{ij}^\beta}\]
Para ello calculamos primero la distancia del Ikea al centro de cada barrio
cc <- coordinates(mapp) dd <- spDistsN1(cbind(cc[,1],cc[,2]),c(lon_Ikea,lat_Ikea),longlat = T)
Hacemos el fit (en logaritmos ya que es una función no lineal)
fit <- lm(log(nusers) ~ log(dd), data=mapp@data) summary(fit)
## ## Call: ## lm(formula = log(nusers) ~ log(dd), data = mapp@data) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.55941 -0.46282 -0.04791 0.46385 1.78361 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.6864 0.4018 11.66 < 2e-16 *** ## log(dd) -1.4617 0.1757 -8.32 3.27e-13 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.7135 on 106 degrees of freedom ## Multiple R-squared: 0.3951, Adjusted R-squared: 0.3893 ## F-statistic: 69.22 on 1 and 106 DF, p-value: 3.274e-13
donde vemos que \(\beta \approx\) 1.46
Pero el numero de personas también depende de la población del barrio, asi que vamos a suponer una ley de la forma
\[N_{ij} \sim \frac{P_i^\alpha}{d_{ij}^\beta}\]
con lo que el modelo sería
fit1 <- lm(log(nusers) ~ log(pob) + log(dd), data=mapp@data) summary(fit1) ## ## Call: ## lm(formula = log(nusers) ~ log(pob) + log(dd), data = mapp@data) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.44146 -0.44397 -0.02294 0.41141 1.81288 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 1.5928 1.1091 1.436 0.15393 ## log(pob) 0.3045 0.1023 2.977 0.00361 ** ## log(dd) -1.4454 0.1696 -8.522 1.24e-13 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.6884 on 105 degrees of freedom ## Multiple R-squared: 0.4421, Adjusted R-squared: 0.4315 ## F-statistic: 41.61 on 2 and 105 DF, p-value: 4.927e-14
Dado que el fit es a un número de usuarios (variable discreta), probamos un GLM poissoniano con link logaritmico
\[\log (\lambda_{ij}) \sim \alpha log(P_i) - \beta log(d_{ij})\]
fit2 <- glm(nusers ~ log(pob) + log(dd),data=mapp@data,family=poisson(link="log")) summary(fit2) ## ## Call: ## glm(formula = nusers ~ log(pob) + log(dd), family = poisson(link = "log"), ## data = mapp@data) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -4.5507 -1.1025 -0.3751 0.5127 5.7009 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -2.73895 0.77576 -3.531 0.000415 *** ## log(pob) 0.79023 0.07114 11.107 < 2e-16 *** ## log(dd) -1.62842 0.06308 -25.816 < 2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 1123.18 on 107 degrees of freedom ## Residual deviance: 296.49 on 105 degrees of freedom ## AIC: 657.59 ## ## Number of Fisher Scoring iterations: 5
Los resultados anteriores nos dan una forma funcional para el número de usuarios del barrio \(i\) que visitan el centro comercial \(j\):
\[N_{ij} \sim \frac{P_i^\alpha}{d_{ij}^\beta}\].
Vamos a utilizarla para determinar las zonas de atracción suponiendo que esa ley vale para todos los centros comerciales
nij_Huff_model <- function(lon1,lat1,lon2,lat2,pob,alpha,beta){
dd <- spDistsN1(cbind(lon1,lat1),c(lon2,lat2),longlat = T)
prob_Huff_model <- pob^alpha/dd^beta
prob_Huff_model
}
centros_sel <- centros[centros$ETIQUETA=="Ikea",]
coord_barrios <- coordinates(mapb)
nij <- data.frame()
for(i in 1:nrow(centros_sel)){
lon1 <- centros_sel$lon[i]
lat1 <- centros_sel$lat[i]
nij <- rbind(nij,nij_Huff_model(coord_barrios[,1],coord_barrios[,2],lon1,lat1,
mapb@data$pob,0.31,1.44))
}
colnames(nij) <- mapb@data$GEO
head(nij[,1:8])
## 079011 079012 079013 079014 079015 079016 079021 ## 1 0.4878643 0.5824864 0.3449798 0.3780751 0.4858795 0.3222756 0.5371665 ## 2 0.3829762 0.4738342 0.3221769 0.4014899 0.4776114 0.2876585 0.3530303 ## 3 0.6881971 1.0098479 0.6549087 0.7006212 0.7604057 0.5327172 0.7002193 ## 079022 ## 1 0.5884984 ## 2 0.4077100 ## 3 0.9414437
Para cada barrio, nos quedamos con el centro que da más número de visitantes. Eso define la zona de atracción. Y la visualizamos
centro_most <- apply(nij,2,which.max) mapb@data$centro_most <- centro_most factpal <- colorFactor(topo.colors(3), as.factor(mapb@data$centro_most)) leaflet(mapb) %>% addTiles() %>% addPolygons(weight=2,color = ~factpal(as.factor(centro_most)))
Repetir el análisis de las zonas de atracción suponiendo que el IKEA de Alcorcón cierra.
Repetir el análisis para los centros de Leroy Merlin